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DEVELOPMENT OF CURVED-PLATE ELEMENTS FOR THE EXACT 
BUCKLING ANALYSIS OF COMPOSITE PLATE ASSEMBLIES INCLUDING 

TRANSVERSE-SHEAR EFFECTS 

ABSTRACT 

The analytical formulation of curved-plate non-linear equilibrium equations including 
transverse-shear-deformation effects is presented. The formulation uses the principle of 
virtual work. A unified set of non-linear strains that contains terms from both physical 
and tensorial strain measures is used. Linearized, perturbed equilibrium equations 
(stability equations) that describe the response of the plate just after buckling occurs are 
then derived after the application of several simplifying assumptions. These equations 
are then modified to allow the reference surface of the plate to be located at a distance \ 
from the centroidal surface. The implementation of the new theory into the VICONOPT 
exact buckling and vibration analysis and optimum design computer program is described 
as well. The terms of the plate stiffness matrix using both classical plate theory (CPT) 
and first-order shear-deformation plate theory (SDPT) are presented. The necessary steps 
to include the effects of in- plane transverse and in-plane shear loads in the in-plane 
stability equations are also outlined. Numerical results are presented using the newly 
implemented capability. Comparisons of results for several example problems with 
different loading states are made. Comparisons of analyses using both physical and 
tensorial strain measures as well as CPT and SDPT are also made. Results comparing the 
computational effort required by the new analysis to that of the analysis currently in the 
VICONOPT program are presented. The effects of including terms related to in-plane 
transverse and in-plane shear loadings in the in-plane stability equations are also 
examined. Finally, results of a design-optimization study of two different cylindrical 
shells subject to uniform axial compression are presented. 



TABLE OF CONTENTS 


ABSTRACT i 

TABLE OF CONTENTS ii 

LIST OF TABLES iv 

LIST OF SYMBOLS vii 

CHAPTER I 1 

INTRODUCTION 1 

1.1 Purpose of Study 1 

1.2 Literature Review 4 

1.3 Scope of Study 10 

CHAPTER II 12 

ANALYTICAL FORMULATION 12 

2. 1 Plate Geometry, Loadings, and Sign Conventions 12 

2.2 Strain- Displacement Relations 13 

2.3 Equilibrium Equations 16 

2.4 Stability Equations 19 

2.5 Stability Equations Transformed to the Plate Reference Surface 23 

2.6 Constitutive Relations 27 

CHAPTER III 31 

IMPLEMENTATION INTO VICONOPT 31 

3. 1 Simplifications to the Theory 31 

3.2 Continuity of Rotations at a Plate Junction 32 

ii 



3.3 Derivation of the Curved-Plate Stiffness Matrix 


36 


3.4 The Wittrick-Williams Eigenvalue Algorithm 42 

CHAPTER IV 44 

NUMERICAL RESULTS 44 

4. 1 Convergence of the Segmented-Plate Approach 44 

4.2 Buckling of Curved Plates With Widely Varying Curvatures 46 

4.3 Buckling of an Unsymmetrically Laminated Curved Plate 48 

4.4 Effect of N 22 Terms in the In-Plane Stability Equations 49 

4.5 Design Optimization of a Cylindrical Shell Subject to Uniaxial Compression .... 53 

CHAPTER V 56 

CONCLUDING REMARKS 56 

REFERENCES 59 

APPENDIX A 63 

MATRICES FOR DETERMINING CHARACTERISTIC ROOTS 63 


iii 



LIST OF TABLES 


Table 1. Material properties for boron/epoxy plies and 2024 aluminum (SI units). 71 

b 2 

Table 2. Geometric parameters used to vary the curvature parameter, . 71 

R 2 t 

Table 3. Critical value of stress resultant N u for buckling of a symmetrically 

laminated long curved plate with clamped longitudinal edges. 72 

Table 4. Critical value of stress resultant N 22 for buckling of a symmetrically 

laminated long curved plate with clamped longitudinal edges. 73 

Table 5. Critical value of stress resultant N 12 for buckling of a symmetrically 

laminated long curved plate with clamped longitudinal edges. 74 

Table 6. Critical value of stress resultant N u = N 22 = N 12 for buckling of a 

symmetrically laminated long curved plate with clamped longitudinal 
edges. 75 

Table 7. Critical value of stress resultant N u for buckling of an unsymmetrically 

laminated long curved plate with simply supported longitudinal edges. 76 

Table 8. Critical value of stress resultant N 12 for buckling of an unsymmetrically 

laminated long curved plate with simply supported longitudinal edges. 76 

Table 9. Material properties for aluminum and Korex tm honeycomb core (English 

engineering units). 77 

Table 10. Critical value of hoop stress resultant N 22 for buckling of a long, isotropic 

cylinder subject to uniform external compression (results are in lbs/in.). 77 

Table 11. Design-optimization results for a honeycomb-sandwich cylinder subjected 

to N n loading. 77 

Table 12. Design-optimization results for a solid-wall cylinder subjected to N n 

loading. 78 


IV 



LIST OF FIGURES 


Figure 1.1. .Typical longitudinally stiffened plate structures. 79 

Figure 1.2 Segmented representation of curved-plate geometry currently used by 

VICONOPT. " 79 

Figure 2. 1 Curved-plate geometry and sign convention for buckling displacements, 

rotations, moments, and forces. 80 

Figure 2.2 Sign convention for applied in-plane loads and relation of reference 

surface to centroidal surface. 8 1 

Figure 2.3 Curved-laminate geometry. 82 

Figure 3.1 Displacements and rotations at a typical plate junction. 83 

Figure 4. 1 Long isotropic (aluminum) cylinder subjected to uniaxial compression. 84 

Figure 4.2 Convergence of VICONOPT segmented-plate results as a function of the 

number of segments used in the approximation. 85 

Figure 4.3 Normalized CPU time requirements for the segmented-plate approach 

as a function of the number of segments used in the approximation. 86 

Figure 4.4 Positive applied in-plane loads on a long curved plate. 87 

Figure 4.5 Symmetrically laminated long curved plate with clamped longitudinal 

edges subjected to applied in-plane loads. 87 

Figure 4.6 Critical value of stress resultant N u for buckling of a symmetrically 

laminated curved plate with clamped longitudinal edges. 88 

Figure 4.7 Critical value of stress resultant N 22 for buckling of a symmetrically 

laminated curved plate with clamped longitudinal edges. 89 

Figure 4.8 Critical value of stress resultant N 12 for buckling of a symmetrically 

laminated curved plate with clamped longitudinal edges. 90 

Figure 4.9 Critical value of stress resultants N u = N 22 = N ]2 for buckling of a 

symmetrically laminated curved plate with clamped longitudinal edges. 91 

Figure 4.10 Unsymmetrically laminated aluminum and boron/epoxy (B/E) curved 

plate with simply supported edges subjected to applied in-plane loads. 92 

Figure 4. 1 1 Critical value of stress resultant N u for buckling of an unsymmetrically 
laminated aluminum and boron/epoxy (B/E) curved plate with simply 
supported longitudinal edges. 93 

Figure 4. 12 Critical value of stress resultant N ]2 for buckling of an unsymmetrically 
laminated aluminum and boron/epoxy (B/E) curved plate with 
supported longitudinal edges. 94 


v 



Figure 4.13 Isotropic (aluminum) long cylindrical tube subjected to uniform external 

pressure loading. 95 


Figure 4. 14 Cylindrical shell subjected to uniform axial compression (N n loading). 96 

Figure 4.15 Optimized cylinder mass as a function of the applied loading for a 

cylindrical shell. 97 


vi 



LIST OF SYMBOLS 


A extensional stiffness matrix 

a upper half of the eigenvectors of matrix R , associated with displacements 

B coupling stiffness matrix 

B, C, E, 

F, G, H coefficients used to select physical or tensor strains 
b lower half of the eigenvectors of matrix R , associated with forces 

b plate width (arc length) 

C matrix whose columns contain the eigenvectors of matrix R 

c single eigenvector of matrix R 

D bending stiffness matrix 

d vector of displacement amplitudes at the two edges of a plate 

ds arc length of a line element in a body before deformation 

ds arc length of a line element in a body after deformation 

E H Young’s modulus in the i-i direction 

E matrix used to define vector d , see Eq. (3.19) 

F matrix used to define vector f, see Eq. (3.20) 

f vector of force amplitudes at the two edges of a plate 

f body forces 

G 12 in-plane shear stiffness 

G 13 , G 23 transverse shear stiffnesses 

hy coefficient of the partially inverted constitutive relations, see Eqs. (3.11) and 

(3.12) 

I identity matrix 

I moment of inertia 

i imaginary number, square root of - 1 

vii 



K 

k 

M n , M 22 , M 12 
m n , m 22 , m 12 
m n ,m 22 ,m 12 

N N N 

1 Ml’ 1 M2’ 1 M2 

n n , n 22 , n 12 

n ll ’ n 22 ’ n l2 
^ 22 ’ °12 
n i 

P 

P 

P 2 ’P 3 

Q 

Q 

QmQ 2 

qi. q 2 ’ 

qi.q2 

q 2 

R 


plate stiffness matrix 

transverse- shear compliance matrix 

applied (prebuckling) moment resultants 

perturbation values of moment resultants just after buckling has occurred 
moment resultants 

applied (prebuckling) stress resultants 

perturbation values of stress resultants just after buckling has occurred 
stress resultants 

effective forces per unit length at an edge § 2 = constant 

number of layers in a general curved laminate 

coefficient matrix of the set of first-order plate differential equations, 

see Eq. (3.14) 

applied external pressure 

perturbation values of the applied pressure load in the buckled state in 
the § 2 - and ^-directions 

lamina reduced stiffness matrix 

lamina reduced transformed stiffness matrix 

applied (prebuckling) shear stress resultants 

perturbation values of shear stress resultants just after buckling has 

occurred 

shear stress resultants 

effective transverse shear force per unit length at an edge § 2 = constant 

matrix whose eigenvalues are the characteristic roots of the plate 
differential equations, see Eq. (3.16b) 
radii of lines of principal curvature 

viii 


R p R 2 



T 


coefficient matrix of the set of first-order plate differential equations, see Eq. 
(3.14) 

surface tractions 


T, 

t plate thickness 

Uj , U 2 prebuckling displacements 

Uj , u 2 perturbation values of displacements just after buckling has occurred 
v volume 

w normal displacement in the ^-direction 

Z vector of the forces and displacements in the plate 

z vector containing the amplitudes of the forces and displacements in the plate 

assuming a sinusoidal variation in the -direction 

z c distance from the plate centroidal surface to the plate reference surface 

z k distance from laminate reference surface to the kth layer in the laminate 

Greek 

a,, a 2 Lame parameters 

(3 angle included by a curved plate 

e vector containing strains e H , e 22 , and y ]2 

e n , e 22 in-plane direct strains 

e 12 , Y 12 in-plane shear strains 

e 13 , Y 13 transverse shear strains 

e 23 , Y23 transverse shear strains 

(Jr, , (j ) 2 rotations 


IX 



<t>, 


K n , K 22 


k 


0 k 


p 

a 




22 


b l2 


l„i 2 .i 


3 


rotation about the normal to the plate middle surface 
middle surface changes in curvatures 
middle surface twisting curvature 
half wavelength of buckling mode 
Poisson’s ratio 

angular orientation of ply k in a laminate with respect to the laminate 

coordinate system 

density 

vector containing stresses o u , o 22 , and x 12 
in-plane direct stresses 
in-plane shear stress 

coordinate measures in the 1-, 2-, and 3-directions, respectively 


Subscripts and Superscripts 

cr critical value for buckling 

k kth layer in a laminated composite plate 

n normal to middle surface 

1, 2, 3 1-, 2-, and 3-directions, respectively 

0 value at centroidal surface 


x 



CHAPTER I 


INTRODUCTION 
1.1 Purpose of Study 

Longitudinally stiffened plate structures occur frequently in aerospace vehicle 
structures. These structures can typically be represented by long, thin, flat or curved 
plates that are rigidly connected along their longitudinal edges, see Figure 1.1. The 
designs for these structures often exploit the increased structural efficiency that can be 
obtained by the use of advanced composite materials. Therefore, the plates used to 
represent the structure may consist of anisotropic laminates. The buckling and vibration 
behavior of this type of structure must be understood to design the structure. 
Additionally, to satisfy the current demands for more cost-effective and structurally 
efficient aerospace vehicles, these structures are frequently optimized to obtain an 
optimal design that satisfies either buckling or vibration constraints or a combination of 
these two constraints. There is a need for analytical tools that can provide the analysis 
capability required to optimize panel designs. 

The VICONOPT computer code [1] is an exact analysis and optimum design 
program that includes the buckling and vibration analyses of prismatic assemblies of flat, 
in- plane- loaded anisotropic plates. The code also includes approximations for curved and 
tapered plates, discrete supports, and transverse stiffeners. Anisotropic composite 
laminates having fully populated A , B and D stiffness matrices may be analyzed. Either 
classical plate theory (CPT) or first-order transverse-shear-deformation plate theory 
(SDPT) may be used [2], The analyses of the plate assemblies assume a sinusoidal 
response along the plate length. The analysis used in the code is referred to as “exact” 
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because it uses stiffness matrices that result from the exact solution to the differential 


equations that describe the behavior of the plates. 

Currently, VICONOPT approximates a curved plate by subdividing it into a series of 
flat-plate segments that are joined along their longitudinal edges to form the complete 
curved-plate structure, see Figure 1.2. This procedure is analogous to the discretization 
approach used in finite element analysis. The code uses exact stiffnesses for the flat-plate 
segments and enforces continuity of displacements and rotations at the segment 
connections. Thus, the analyst must ensure that an adequate number of flat-plate 
segments is used in the analysis. The next logical step in the development of the 
VICONOPT code is to eliminate the need to approximate curved-plate geometries by 
flat-plate segments by adding the capability to analyze curved-plate segments exactly. 
By adding this capability, the accuracy of the solutions can be improved. Furthermore, 
since the curvature of a plate is modeled directly, there will be no need to determine if a 
sufficient amount of flat-plate segments have been used to model the curved plate. 
Another benefit of adding this capability is that the computational efficiency of the code 
will be improved since only one stiffness calculation for the entire curved plate is 
required, rather than the several that are currently required for the individual flat plates 
that are used to approximate the curved plate. This improvement in computational 
efficiency is important for structural optimization. In this report, the capability to analyze 
curved-plate segments exactly has been added to the VICONOPT code. The present 
report will describe the methodology used to accomplish this enhancement of the code 
and will present results obtained utilizing this new capability. 

The procedure used in the present report is an extension of the procedure described in 
[2], This procedure involves deriving the appropriate differential equations of 
equilibrium for the analysis of fully anisotropic curved plates, including transverse-shear- 
deformation effects. These coupled equations are of eighth-order if transverse-shear 
effects are neglected, and of tenth-order if transverse-shear effects are included. For the 
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analysis of flat plates, the coupling of these equations occurs through the laminate 
extension-bending B matrix; however, coupling can also be produced by including 
curvature terms in the equilibrium equations. The numerical solution technique that was 
developed in [2] to solve such systems of equations will apply for either type of coupling, 
and the stiffnesses of the plates are derived from the numerical solution to these 
equations. 

Several features have been added to the VICONOPT code as part of the present 
report. The current version of VICONOPT only analyzes flat-plate elements based on a 
tensorial strain-displacement relation. However, the choice of strain-displacement 
relations can affect the contribution of prebuckling forces in curved plates. Therefore, a 
unified set of nonlinear strain-displacement relations that contains terms from both 
physical and tensorial strain measures is used to derive the plate equilibrium equations. 
The unified set of strains is used throughout the derivation of the equilibrium equations, 
and the selection of either physical or tensorial strains is achieved by appropriately 
setting coefficients in the equilibrium equations equal to one or zero. The option to use 
physical strain-displacement relations for the analysis of flat plates is included as well. 
Another addition is the treatment of the effects of in-plane transverse and in-plane shear 
loadings in the in-plane equilibrium equations. These effects are currently ignored in the 
VICONOPT code (see [1]). In the present report, an in-plane transverse loading, denoted 
N 22 , is a loading that acts perpendicular to the longitudinal edges of the plate. The 
present study has added the option to include the effects of these loadings in the in-plane 
equilibrium equations. Finally, either CPT or SDPT may be used. The SDPT used in 
VICONOPT and in the present report uses the usual first-order assumption that straight 
lines originally normal to the centroidal surface are assumed to remain straight and 
inextensional but not necessarily normal to the centroidal surface during deformation of 
the plate. All of these features have been implemented such that they are available for 
use in the analysis of both flat and curved plates. 
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1.2 Literature Review 


The buckling and vibration analysis of assemblies of prismatic plates has received a 
great deal of attention over the last thirty years. One method of analysis for this class of 
structure that has been studied extensively is the finite-strip method, FSM [3]. A popular 
application of this method involves determining a stiffness matrix for each individual 
plate in the assembly and then assembling those individual matrices into a global stiffness 
matrix for use in determining the response of the entire structure. This method is 
therefore analogous in form to the finite element method [4], The main difference 
between the two methods is that the finite element method discretizes the individual 
plates into elements in both the longitudinal and transverse directions. The stiffness 
matrix for each individual element is then calculated and assembled into a global stiffness 
matrix. In the FSM, the response of the plate in the longitudinal direction is represented 
as a continuously differentiable smooth series that satisfies the boundary conditions at the 
two ends of the plate. Therefore, discretization of the structure is only required to be 
performed in the transverse direction, and depending on the method being used, 
discretization of the individual plates may or may not be required [3]. 

The work in the area of finite strip analysis of assemblies of prismatic plates may be 
broadly classified based upon different characteristics of the analysis method used. One 
classification distinguishes whether the properties of the individual plates are derived by 
direct solution to the equations of equilibrium or by application of potential energy or 
virtual work principles, i.e., exact versus approximate methods. Another classification 
distinguishes whether classical plate theory (CPT) or first-order shear-deformation plate 
theory (SDPT) is used in the analysis. Finally, a distinction may be made as to whether 
or not complex quantities are used in the development of the individual stiffness matrices. 
A review of the literature in the area of finite strip analysis methods is presented below. 
Approximate methods are discussed separately from exact methods. 
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The approximate FSM was first proposed for the static analysis of plate bending by 
Cheung in 1968 [5], The approximate FSM involves subdividing each plate into a series 
of finite-width strips that are linked together at their longitudinal edges in a manner 
similar to that depicted in Figure 1.2. Separate expressions for in-plane and out-of-plane 
displacements as well as rotations about the in-plane x and y axes over the middle surface 
of each strip are assumed. Each of these fundamental quantities are expressed as a 
summation of the products of longitudinal series and transverse polynomials [3]. The 
longitudinal series are typically sinusoidal and are selected to satisfy displacement 
conditions at the transverse edges of each strip that match the desired plate boundary 
conditions along those edges. The potential energy of an individual finite strip is then 
evaluated, and the total potential energy of the plate is obtained by summing the potential 
energies of the individual strips. Following the application of any appropriate zero- 
displacement boundary conditions at the longitudinal edges, the potential energy is 
minimized with respect to each plate degree of freedom to generate the equilibrium 
equations for the plate. Displacements are then calculated for a given loading condition 
using this system of equations. 

The analysis of [5] utilized CPT for the static bending analysis of isotropic plates. In 
1971, Cheung and Cheung [6] applied the approximate FSM to the analysis of natural 
vibrations of thin, flat-walled structures with different combinations of the standard edge 
boundary conditions (i.e., clamped, simply supported, or free). Their analysis was based 
upon CPT and the displacements in the longitudinal direction were approximated using 
the normal modes of Timoshenko beam theory to allow for various boundary conditions 
on the transverse edges. 

Przemieniecki [7] used an approximate FSM based upon CPT to calculate the initial 
buckling of assemblies of flat plates subjected to a biaxial stress state. This method only 
considered local buckling modes since it assumed that the line junctions between plates 
remained straight during buckling. Plank and Wittrick [8] extended the work of 
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Przemieniecki by considering global as well as local modes and by admitting a more 
general loading state that included uniform transverse and longitudinal shear stress and 
longitudinal direct stress that varies linearly across the width of the plate. When in-plane 
shear loading is present, a spatial phase difference occurs between the perturbation forces 
and displacements which occur at the edges of the plates during buckling. This phase 
difference causes skewing of the nodal lines and is accounted for in [8 ] by defining the 
magnitude of these quantities using complex quantities. This method is referred to as a 
complex finite strip method. 

In 1977, Dawe [9 and 10] used an approximate FSM based upon CPT for the static 
and linear buckling analysis of curved-plate assemblies. The plates studied were 
isotropic, and in-plane shear loads were not allowed. Morris and Dawe extended this 
analysis to study the free vibration of curved-plate assemblies in 1980 1 1 1]. 

All of the analyses discussed thus far have been based upon CPT. In 1978, Dawe 
[12] presented an approximate FSM based upon SDPT [13] for the vibration of isotropic 
plates with a pair of opposite edges simply supported. Roufaeil and Dawe [14] and 
Dawe and Roufaeil [15] extended this analysis to the vibration and buckling, 
respectively, of isotropic and transversely isotropic plates with general boundary 
conditions. The latter two analyses admitted the general boundary conditions through the 
use of the normal modes of Timoshenko beam theory, as was done in [6]. 

In 1986, Craig and Dawe [16] considered the vibration of single symmetrically 
laminated plates using an approximate FSM based upon SDPT. Dawe and Craig [17] 
then extended this analysis to study the buckling of single symmetrically laminated plates 
subject to uniform shear stress and direct in-plane stress. This analysis allowed for 
anisotropic material properties. General boundary conditions were once again admitted 
through the use of the normal modes of Timoshenko beam theory. The analysis of [17] 
was extended in 1987 to the vibration of complete plate assemblies [18]. However, it 
was shown in this work that the problem size increased dramatically as attempts to 
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increase the accuracy of the solution were made by further subdivision of the component 
plates. 

In 1988, Dawe and Craig [19] presented a complex FSM based upon SDPT for the 
buckling and vibration of prismatic plate structures in which the component plates could 
consist of anisotropic laminates and could be subject to in-plane shear loads. This work 
also made use of substructuring to create “superstrips” that eliminated the internal 
degrees-of-freedom from each component plate. This analysis was later extended to 
consider finite-length structures [20 and 21] and to add multi-level substructuring to 
couple several “superstrips” to further decrease the problem size. Dawe and Peshkam 
[22] also developed a complementary analysis to that presented in [20 and 21] for long 
plate structures. Analyses using both SDPT and CPT were presented. This work also 
added the capability to define eccentric connections of component plates. 

Wittrick laid the groundwork for the exact FSM in 1968 [23]. The basic assumption 
in this work is that the deformation of any component plate varies sinusoidally in the 
longitudinal direction. Using this assumption, a stiffness matrix may be derived that 
relates the amplitudes of the edge forces and moments to the corresponding edge 
displacements and rotations for a single component plate. For the exact FSM, this 
stiffness matrix is derived directly from the equations of equilibrium that describe the 
behavior of the plate. In [23], Wittrick developed an exact stiffness matrix for a single 
isotropic, long flat plate subject to uniform axial compression. His analysis used CPT. 
Wittrick and Curzon [24] extended this analysis to account for the spatial phase 
difference between the perturbation forces and displacements which occur at the edges of 
the plate during buckling due to the presence of in-plane shear loading. This phase 
difference is accounted for by defining the magnitude of these quantities using complex 
quantities. Wittrick [25] then extended his analysis to consider flat isotropic plates under 
any general state of stress that remains uniform in the longitudinal direction (i.e., 
combinations of bi-axial direct stress and in-plane shear). A method very similar to that 
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described in [23] was presented by Smith in 1968 [26] for the bending, buckling, and 
vibration of plate-beam structures. 

In 1972, Williams [27] presented two computer programs, GASVIP and VIPAL to 
compute the natural frequencies and initial buckling stress of prismatic plate assemblies 
subjected to uniform longitudinal stress or uniform longitudinal compression, 
respectively. GASVIP was used to set up the overall stiffness matrix for the structure, 
and VIPAL demonstrated the use of substructuring. In 1974, Wittrick and Williams [28] 
first reported on the VIPASA computer code for the buckling and vibration analyses of 
prismatic plate assemblies. This code allowed for isotropic or anisotropic plates as well 
as a general state of stress (including in-plane shear). The complex stiffnesses described 
in [8] were incorporated, as well as allowances for eccentric connections between 
component plates. This code also incorporated an algorithm, referred to as the Wittrick- 
Williams algorithm, for determining any natural frequency or buckling load for any given 
wavelength [29]. The development of this algorithm was necessary because the complex 
stiffnesses described above are transcendental functions of the load factor and half 
wavelength of the buckling modes of the structure. The eigenvalue problem for 
determining natural frequencies and buckling load factors is therefore transcendental. 
Further discussion of the Wittrick-Williams algorithm will be presented in Chapter III. 

In 1973, Viswanathan and Tamekuni [30 and 31] presented an exact FSM based upon 
CPT for the elastic stability analysis of composite stiffened structures subjected to biaxial 
inplane loads. The structure is idealized as an assemblage of laminated plate elements 
(flat or curved) and beam elements. The analysis assumes that the component plates are 
orthotropic. The transverse edges are assumed to be simply supported, and any 
combination of boundary conditions may be applied to the longitudinal edges. The 
analysis was included in an associated computer code, BUCLAP2. Viswananthan, 
Tamekuni, and Baker extended this analysis in [32] to consider long curved plates 
subject to any general state of stress, including in-plane shear loads. Anisotropic material 
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properties were also allowed. This analysis utilized complex stiffnesses as described in 
[8]. The analyses described in [26, 28, and 32] are very similar. The differences between 
the three are discussed in [28]. 

When applied in-plane shear loads or anisotropy is present, the assumption of a 
sinusoidal variation of deformation in the longitudinal direction is only exact for 
structures that are infinitely long. Significant errors for structures of finite length can 
occur due to the skewing of nodal lines. In 1983, Williams and Anderson [33] presented 
modifications to the eigenvalue algorithm described in [29]. The modifications presented 
in [33] allowed the buckling mode corresponding to a general loading to be represented 
as a series of sinusoidal modes in combination with Lagrangian multipliers to apply point 
constraints at any location on those edges. Each sinusoidal mode is represented by an 
exact stiffness matrix. This technique allows infinitely long structures supported at 
repeating intervals with anisotropy or applied in-plane shear loads to be analyzed. Thus, 
a panel supported at its transverse edges is approximated by one with a series of point 
supports along those edges. These modifications formed the basis for the computer code 
VICON (VIpasa with CONstraints) described in [34], However, the analysis capability 
of VICON was limited to plates analyzed with CPT having a zero B matrix. The VICON 
code was later modified to include structures supported by Winkler foundations [35]. An 
optimum design feature was also added in 1990 [36 and 37], and the VICONOPT 
(VICON with OPTimization) code was introduced. 

Anderson and Kennedy [2] incorporated SDPT into VICONOPT in 1993. A 
numerical approach to obtain exact plate stiffnesses that include the effects of transverse- 
shear deformation was presented. The generality of VICONOPT was also expanded in 
[2] to allow for the analysis of laminates with fully populated A , B , and D stiffness 
matrices. 
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1.3 Scope of Study 

The analytical formulation of the curved-plate non-linear equilibrium equations 
including transverse-shear-deformation effects are presented in Chapter II. A unified set 
of non-linear strains that contains terms from both physical and tensorial strain measures 
is used. The equilibrium equations are derived using the principle of virtual work 
following the method presented by Sanders [38 and 39]. Linearized, perturbed 
equilibrium equations that describe the response of the plate just after buckling occurs are 
then derived after the application of several simplifying assumptions. Modifications to 
these equations that allow the reference surface of the plate to be located at a distance \ 
from the centroidal surface are then made. 

In Chapter III, the implementation of the new theory into the VICONOPT code is 
described. A derivation of the terms of the plate stiffness matrix using MATHEMATICA 
[40] is presented. The form of these terms for both CPT and SDPT is discussed. The 
necessary steps to include the effects of in-plane transverse and in-plane shear loads in 
the in-plane equilibrium equations are also outlined. 

In Chapter IV, numerical results are presented using the newly implemented 
capability. A convergence study using the current segmented-plate approach in 
VICONOPT is performed for a simple example problem to obtain baseline results for use 
in future comparisons. Results comparing the computational effort required by the new 
analysis to that of the analysis currently in the VICONOPT program are also presented. 
Comparisons of results for several example problems with different loading states are 
then made. Comparisons of analyses using both physical and tensorial strain measures as 
well as CPT and SDPT are made. The effects of including terms related to in-plane 
transverse and in-plane shear loads in the in-plane stability equations are also examined. 

In Chapter V, the characteristics of the newly implemented curved-plate elements in 
VICONOPT is presented. A brief summary of the effects of several analytical features 
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that have been implemented into VICONOPT is given. Finally, potential future work in 
this area is discussed. 
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CHAPTER II 


ANALYTICAL FORMULATION 

In this chapter, the non-linear equilibrium equations are derived for a curved plate 
including transverse-shear effects. A unified set of non-linear strains that contain terms 
from both physical and tensorial strain measures is used. The equilibrium equations are 
derived using the principle of virtual work following the method presented by Sanders 
[38 and 39]. Linearized stability equations that describe the response of the plate just 
after buckling occurs are then derived following the application of several simplifying 
assumptions. Modifications to these equations that allow the reference surface of the 
plate to be located at a distance z c from the centroidal surface are then made. 

2.1 Plate Geometry, Loadings, and Sign Conventions 

The geometry of the basic plate element being studied is given in Figure 2.1. This 
figure depicts the orthogonal curvilinear coordinate system (§,, § 2 , § 3 ) used in the present 

analysis. The |j- and § 2 -axes shown in the figure are along lines of principal curvature 

and they have radii of curvature R, and R 2 , respectively. The § 2 -axis is normal to the 

middle surface of the plate. The first fundamental form of the plate middle surface is 
given by 

ds 2 =a 2 d§ 2 +a 2 d ^2 (2.1) 

where a, and a 2 are the Lame parameters. The coordinates and % are measured as arc 
lengths along the and | 2 -axes, respectively. The result of measuring the coordinates 
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in this manner is that dj = a 2 , = 1. The sign conventions for buckling displacements, 

moments, rotations, and forces are also shown in Figure 2.1. The sign convention for the 
applied in-plane loadings being considered and the relation of the reference surface of the 
plate to the centroidal surface of the plate are shown in Figure 2.2. Note that that 
centroidal surface can be offset from the reference surface by a distance z c . The 
centroidal surface is defined to be located at the centroid of the face of the panel that is 

normal to the ^,-axis. The loading N 22 shown in this figure is referred to in the present 

report as an in-plane transverse loading. 

2.2 Strain-Displacement Relations 

The nonlinear strain-displacement relations used for the present study are given by 
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(2.2e) 


U-) 

2e 23 = Y 23 = W ’2 - “ -( t ) 2 


where the following notation for partial derivatives is used: -^- = Ujj. The 

displacement quantities in Eqs. (2.2a) through (2.2e) are displacements of the centroidal 
surface of the plate. The constants B, C, E, F, and H are set equal to one and G is set 
equal to zero in Eqs. (2.2a) through (2.2e) to use tensorial strain measures. The constants 
B, E, and G are set equal to one and C, F, and H are set equal to zero to use physical 
strain measures. Note that the linear portions of the tensorial and physical strain 
measures are identical. To obtain Donnell theory from the strain-displacement relations 
in Eqs. (2.2a) through (2.2e) the constants B, C, E, F, G_, and H must be set equal to zero, 


U, Uo 

and all terms involving the quantities — and must be neglected. Sander’s theory 

Rj R 2 

[39] may be obtained by setting the constants B, C, E, F, G_, and H equal to zero and 


1 2 

adding the term — <|) n to Eqs. (2.2a) and (2.2b), where <j) n is the rotation about the normal 


to the plate middle surface. 

The tensorial strain measures used in the present study are those of Novozhilov [41], 
These strains are obtained by taking the difference between the square of the arc length of 
a line element in a body after deformation, (ds*) 2 , and before deformation (ds) 2 . The 

tensorial strain measures, e jk , are defined by the relationship 



Ejkd§jd§ k 


i, j = 1, 3 


(2.3) 


The repeated indices in Eq. (2.3) indicate summation over i and j. The physical strain 
measures are strains that can be measured in the laboratory. The physical strains used in 
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the present report are derived in a manner similar to that presented by Stein in [42] and 
they were communicated to the author in lines of curvature coordinates by Dr. Michael P. 
Nemeth. 1 . Physical extensional strains are defined as the ratio of the change in arc length 
of a line element in a body, ds\ to the original length of that line element, ds, 

(ds*) -(ds)j 

e ;; = j = 1,2 (no summation) (2.4a) 

(ds)j 

Physical shearing strains are defined as the change in the angles between three line 
elements that are orthogonal before deformation and are oriented in the direction of three 

A 

unit vectors, ej, after deformation. The physical shearing strains are defined by the 
following expressions 


si n (y 12 ) ** Y 12 = (2-4b) 

sin (Yj 3 ) 5S Yj3 =ej # e 3 j - 1,2 (2.4c) 


The definitions for the changes in curvatures of the centroidal surface used for both 
theories are 


K 11 = “4>l,l 

(2.5a) 

k 22 = — 4*2,2 

(2.5b) 

K 12 = “( 4 * 1,2 + 4*2,1 ) 

(2.5c) 


These changes in curvatures are equivalent to those given by Sanders in [39] with the 
terms involving rotations about the normal neglected. 


1 Mechanics and Durability Branch, Structures and Materials, NASA Langley Research Center, Hampton, VA, 
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2.3 Equilibrium Equations 

The nonlinear equilibrium equations for the curved plate illustrated in Figures 2. 1 and 
2.2 are derived using the principle of virtual work [43]. This principle states that, if a 
structure in equilibrium is subject to a virtual distortion while remaining in equilibrium, 
then the external virtual work done by the external forces on the structure is equal to the 
internal virtual work done by the internal stresses. The principle of virtual work can 
therefore be written in the form 


jT;6ujds+ /fj6ujdv= JOjjbejjdv (2.6) 

surface volume volume 


The present derivation uses the principle of virtual work in the manner of Sanders [38] 
written in the following form 


ff 

area 


d§ldi; 
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f[ N ll 8u l +Ni 2 Su 2 +Q 1 6w-M 1] 6(|) ] -M 12 6(|)2]d§2 


(2.7) 


- f[Nl26ui + N226U2 + Q26W - M]2&<j>i - M 2 2&<t>2]d?i = 0 


The terms h 12 and m ] 2 are effective stress measures as defined by Sanders in [38]. The 
terms q | and q 2 are also effective stress measures as defined by Cohen in [44]. The 
uppercase terms in Eq. (2.7) are applied loadings on the boundary of the plate. 

Substituting Eqs. (2.2a) through (2.2e) and Eqs. (2.5a) through (2.5c) into Eq. (2.7) 
and integrating by parts results in 
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+ *11,1 + A 12,2 “Ql *t>j + *12,1 + *22,2 ~Al 
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For arbitrary displacements u,, u 2 , w, <j>,, and (|) 2 , the coefficients of the displacements in 

the area integral in Eq. (2.8) are the five equilibrium equations. The coefficients of the 
displacement variables in the first line integral in Eq. (2.8) are the natural boundary 
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conditions for an edge = constant, and the coefficients of the displacement variables in 


the second line integral are the natural boundary conditions for an edge § 2 = constant. 

2.4 Stability Equations 

A set of perturbation equilibrium equations that govern the stability of the plate, 
referred to herein as the stability equations, may now be written by taking the difference 
between the equilibrium equations evaluated for an equilibrium state just prior to 
buckling and an adjacent (perturbed) equilibrium state just after buckling has occurred. 
Let the prebuckling state be represented by: 

5 11 = -Nil* n 22 = “ N 22> n 12 = -N 12> m ll = -M 11> 

A 22 = ~^ 22 i *12 = Ql = “Ql> <12 = “Q2> (2-9) 

Uj, u 2 , w 

The minus signs in the loading terms reflect the sign convention used in which the 
applied loads are opposite in direction to the loads that develop after buckling. Let the 
perturbed state just after buckling has occurred be represented by: 

fi ll = n ll -N 11’ S 22 = n 22 -N 22’ fi 12 = n 12 “ N 12 » 

*11 = m n -M,i, *22 = m 22 - M 2 2 , *12 = m 12 “ M 12 . (2.10) 

qi = qi -Qi, q 2 = q 2 - Q 2 ’ u i + Ui » u 2 + U 2 > w + w 

where the lower case variables are perturbation variables. Taking the difference between 
the two equilibrium states represented by Eqs. (2.9) and (2.10), linearizing the resulting 
equations for the perturbation variables, and applying the following simplifying 
assumptions: 

1) Prebuckling deformations, moments, and transverse- shear stresses are 
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negligible 


2) The in-plane prebuckling stress state is uniform 


yields the following stability equations: 
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The boundary conditions for an edge = constant are 


(2.11a) 


(2.11b) 


(2.11c) 


(2. lid) 
(2. lie) 
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As will be discussed in Chapter III, a sinusoidal variation of displacements and forces is 
assumed in the |j direction. Therefore, these boundary conditions are ignored herein. 

The boundary conditions for an edge fc; 2 = constant are 
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or 


(2.13c) 
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where the terms with a caret ( A ) are effective force quantities per unit length at an edge 
= constant. The effective forces, n 12 ,h 22 > anc ' q 2 are equal to forces in the original 

(undeformed) |j-, ij 2 -, and ^-directions along the longitudinal edges of the plate 

(§ 2 =constant). Introduction of these force quantities facilitates the derivation of the 

stiffness matrix in Chapter III which relates the forces along the longitudinal edges of the 
plate in the original coordinate directions to the corresponding displacements along those 
edges. 

The first three stability equations given in Eqs. (2.11a) through (2.11c) are now 
written in a simplified form using the definitions of the effective forces per unit length 
given in Eqs. (2.13a) through (2.13c) 
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(2.14c) 
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This form of these stability equations will be used herein. Note that Eq. (2.14b) contains 
the perturbation variables n 12 and q 2 . These variables are related to the effective forces, 
n 12 and q 2 , through Eqs. (2.13a) and (2.13c). 

2.5 Stability Equations Transformed to the Plate Reference Surface 

The stability equations given in Eqs. (2.11a) through (2. lie) describe the response at 
the centroidal surface of the plate. A superscript 0 may be added to the displacement 
quantities in these equations to indicate that they are centroidal quantities. These 
equations are now written such that they describe the response at the reference surface of 
the plate, which can be located a distance z c from the centroidal surface, Figure 22. To 
write the stability equations at the reference surface, the following information is used: 

1) The relations of the displacements at the centroidal surface, u° and u°, to the 
displacements at the reference surface, u, and u 2 are: 

uf = U! -z c ty\ (2.15a) 

u 2 =u 2- z c < 1 ) 2 (2.15b) 

2) The relations of the moments at the centroidal surface, m , m ° 2 * and m ° 2 , to 
the displacements at the reference surface, m^ , m 2 2 , and m 12 are: 
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(2.15c) 


m ll = mu - z c n ll 

m 22 = m 22 “ z c n 22 (2.15d) 

m i 2 = m i2 " z c n i2 (2.15e) 

3) The following quantities do not vary with z: 

N ll> N 22 5 N 12 » nn, n 22’ n 12’ Q2’ an d w 

4) The applied in-plane stresses, N„, N 22 , and N 12 act at the centroidal surface. 

Substitution of Eqs. (2.15a) through (2.15e) into Eqs. (2.14a) through (2.14c) and Eqs. 
(2. 1 Id) and (2. 1 le) yields the following equations 
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m ll,l + m 12,2 - z c ( n n,l +n 12 , 2 )-qi =° (2.16d) 

m 12,l +m 22,2 - z c ( n 124 +n 22,2)-q2 =° (2.16e) 

The natural boundary conditions are also rewritten after substitution of Eqs. (2.15a) 
through (2.15e) into Eqs. (2.13a) through (2.13e). For an edge '% 2 = constant, the natural 
boundary conditions become 
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m 22 “ z c n 22 = 0 


(2.17e) 


The last two stability equations, Eqs. (2.16d) and (2.16e), are now rewritten by 
substituting expressions for the quantities (njy + n 122 ) aQ d ( n 12,1 + 11 22,2 ) d ia t can 

be obtained using Eqs. (2.16a) and (2.16b), respectively, and the definitions for the 
effective forces per unit length, Eqs. (2.17a) through (2.17c). The definitions for the 
effective forces are needed since the terms n 12 and n 22 that appear in the two above are the 
perturbation values, not the effective forces. Substitution of the expressions for the two 
quantities above into Eqs. (2.16d) and (2.16e), respectively, yields the final form of the 
last two stability equations 
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The stability equations in the form given in Eqs. (2.16a) through (2.16c) and Eqs. (2.18a) 
and (2.18b) are those implemented into the VICONOPT code. 

2.6 Constitutive Relations 

The present analysis allows for generally laminated composite materials. The 
geometry of a general, curved laminate is given in Figure 2.3. As shown in the figure, 
the number of layers in the laminate is n,, and the width of the laminate is b. The radius 

of curvature of the § 2 -axis, R 2 is shown in the figure as well. . The radius of curvature of 
the i^-axis, Rj is not shown; however, its direction may be inferred from that of R 2 . The 
lamina coordinate system is the (£,,, ^ 3 ) system and the laminate coordinate system is 

the (|j, | 2 , ^ 3 ) system. The lamina coordinate system is aligned with the principal 

material direction of the lamina, and the laminate coordinate system is aligned with the 
principal geometric directions of the laminate. The coordinate system for the kth lamina 

is oriented at an angle 0 k with respect to the laminate coordinate system. The stress-strain 

relations in the lamina coordinate system for a lamina of orthotropic material in a state of 
plane stress are 
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where the [Q] matrix is referred to as the reduced stiffness matrix for the lamina and is 
defined in [45] in terms of the elastic engineering constants of the lamina. These 
relations may be written in the laminate coordinate system by use of transformation 
matrices as defined in [45]. The transformed relations are 
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(2.20) 
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where the 



matrix is the reduced transformed stiffness matrix for the lamina. Both of 


Eqs. (2.19) and (2.20) may be thought of as stress-strain relations for the kth lamina in a 


multi-layer laminate. Therefore, Eq. (2.20) may be written as 


Mk-[o]kh}k < 2 - 21 > 

The constitutive relations for a thin, elastic laminated composite shell may now be 
defined as 
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where the resultant forces and moments acting on the laminate, {N} and {M}, 
respectively, are defined as 


Nn 

n l 

z k 

a iT 



N 22 

-= 1 

/ < 

°22 

-d§ 3 

(2.23) 

N 12 

k=l 

z k-l 

X 12 , 



Mu' 

n l 

z k 

Oil 



M 22 

= 2 

f ■ 

°22 

- d§ 3 

(2.24) 

M 12 

k=l 

z k-l 

x 12 




28 



where n, is the total number of layers in the laminate. The extensional, coupling, and 
bending stiffness matrices, A, B, and D, respectively, are defined as 

(A, B, D). 2 J [Q] k (l. S 3 . S 3 ) <£3 ( 2 . 25 ) 

k=1 Zk-I 

The analysis in VICONOPT allows for laminates with fully populated A , B , and D 
matrices. 

The constitutive relations for transverse shear used in VICONOPT are those 
presented by Cohen in [44], The constitutive relations for transverse shear are written in 
inverted form as 


J Yl 3 l-f kn 

[Y23j [ k 12 

where |k] is a symmetric 2-by-2 transverse shear compliance matrix whose terms are 
defined in [44j. The terms of the [k] matrix were derived for general, anisotropic, multi- 
layered composite shells and they are a generalization of results for a shell with a 
homogeneous wall for which the transverse shear correction factor for the shear stiffness 
is 5/6. The procedure used in [44] for obtaining the terms of the |k] matrix follows. 
Statically correct expressions of in-plane stresses and transverse-shear stresses were 
derived in terms of the transverse- shear stress resultants and arbitrary constants that were 
interpreted by Cohen as redundant “forces”. The expressions for in-plane stresses were 
obtained using the constitutive relations given in Eq. (2.22) and linear distribution of in- 
plane strains through the wall thickness. The expressions of transverse-shear stresses 

were obtained by integrating in the ^-direction the three-dimensional equilibrium 

equations. The transverse- shear stress resultants were then used to derive an expression 
of the volumetric density of the transverse-shear strain energy. A statically correct 
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expression of the area density of the transverse-shear strain energy was then obtained by 
integrating in the ^-direction this volumetric density. The transverse-shear constitutive 

relations given in Eq. (2.26) were then derived by applying Castigliano’s theorem of least 
work [46] by minimizing the area density of the transverse-shear strain energy with 
respect to the redundant forces mentioned previously. 
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CHAPTER III 


IMPLEMENTATION INTO VICONOPT 


In this chapter, the implementation of the present theory into the VICONOPT code is 
described. Additional simplifications made to the theory are described first. A 

discussion of the use of the transverse-shear strain, y 13 , as a fundamental displacement 

variable in the problem to maintain continuity of rotations at plate junctions is then 
presented. The derivation of an expression for the curved-plate stiffness matrix is 
described. The terms of matrices that are needed to calculate this stiffness matrix were 
obtained using MATHEMATIC A [40], and they are presented in Appendix A. The terms 
for both CPT and SDPT are presented, and the terms that result from the inclusion of 
direct in-plane transverse and in-plane shear loads in the in-plane stability equations are 
specified. As stated previously, the implementation of the curved-plate theory into 
VICONOPT follows very closely the method presented in Reference [2], Therefore, the 
following discussion is necessarily similar to that presented in that reference. 

3.1 Simplifications to the Theory 

Before proceeding with the derivation of the curved-plate stiffness matrix, a 
discussion of several simplifications to be implemented is presented. First, the theory 
implemented into the VICONOPT code considers structures that are prismatic in the 
longitudinal direction. Therefore, for the curved plates being considered in the present 
report, the radius of curvature in the longitudinal direction, R,, is infinite; and any terms 

involving the quantity — — are zero. Although these terms are set equal to zero for the 

R 1 

calculation of the terms of the stiffness matrix, they are retained for completeness in the 
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theory presented in this chapter. Another simplification to the theory involves limiting 
the capability to locate the reference surface a distance z c from the centroidal surface. 
This capability has only been implemented for the case where the effects of N 22 and N 12 
loads in the in-plane stability equations are neglected. The expressions for the stiffness 
terms that result when N 22 and N 12 are included in the in-plane stability equations and z c 
is non-zero are prohibitively long. Therefore, in the derivation to follow, only the 
following two cases are presented: 

1) N 22 and N 12 are included in the in-plane stability equations and z c is zero (i.e., 
reference surface is coincident with the centroidal surface); and , 

2) N 22 and N 12 are neglected in the in-plane stability equations and \ is non-zero (i.e., 
reference surface may be translated from the centroidal surface). 

3.2 Continuity of Rotations at a Plate Junction 

One important issue to be addressed in the analysis of plate assemblies is the 
continuity of rotations at a plate junction. The original VIPASA code is based upon CPT, 
and the theory only treats four degrees of freedom (DOF) at a longitudinal plate edge. 
These DOF are the three displacement quantities, u,, u 2 , and w, and a rotation about the 

^-axis, <j) 2 . Maintaining continuity of these DOF at a typical plate junction is very 

straightforward. However, when SDPT is considered, there are five DOF at a 
longitudinal plate edge. These DOF are the four from CPT as well as an additional 

rotation, <j)j, that results from the inclusion of transverse-shear deformation. Another 
problem that must be addressed is that when two plates are joined together such that one 
is rotated at an arbitrary angle, 0, to the other, rotations about the normals to the 
centroidal surfaces of the two plates must be included to satisfy continuity of rotations. 
This rotation, <j) n , is not accounted for in the present plate theory. The procedure used in 
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VICONOPT to maintain continuity of rotations follows that used by Cohen in [47], This 
procedure introduces the shear strain, y 13 , as a fundamental displacement variable instead 

of the rotation, <(>,. The justification for using this approach is described subsequently. 

The displacements and rotations at a typical plate junction are shown in Figure 3.1. 
The two plates, numbered 1 and 2, are shown viewed along the 1-axis, and it is obvious 
that the Uj displacements are easily matched regardless of the orientation of plate 2. The 

displacements and rotations for which continuity must be maintained are u 2 , w, <|),, and (j) n . 

Upon inspection of Figure 3.1(a), the following expressions for coplanar plates (0 = 0) 
may be written as 


<N (N 

II 

(3.1a) 

= w 2 

(3.1b) 

-♦l 

(3.1c) 

11 

-©- 
3 N> 

(3. Id) 


where the superscripts 1 and 2 refer to the plate numbers. Similarly, upon inspection of 
Figure 3.1(b), the following expressions for 0 = +90° may be written as 


u^ = w 2 (3.2a) 

w 1 = -u 2 (3.2b) 


<t>! = -<t>n (3.2c) 

<l>n = <t>r (3 -2d) 
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Finally, upon inspection of Figure 3.1(c), the following expressions for arbitrary 0 may 


be written as 


2 2 
= u 2 cos0 + w sin0 

(3.3a) 

2 2 
= w cos0 - u 2 sin0 

(3.3b) 

2 2 

= (jij cos0 - <[) n sin0 

(3.3c) 

2 2 

= <|) n cos0 + 4>j sin0 

(3.3d) 


The rotation about the normal of a line element originally directed along the^-axis is 
shown in [48] to be 


du 2 

3 ll 


Using this definition, Eqs. (3.3c) and (3.3d) are written as 


<t>J = <\>] cos0 - U 2 ,i sin0 

12 2 
u 2,l = u 2,l cos0 + <j)j sin0 


(3.4) 


(3.5a) 

(3.5b) 


Using Eqs. (3.3a) and (3.3b) and the definition for y 13 , Eq. (2. 2d), the previous two 
equations may be written as 


y j 3 = cos0 y (3.6a) 

0 = -sin0 y 12 (3.6b) 
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The results shown in Eqs. (3.6a) and (3.6b) indicate that for plates that are not 
coplanar (i.e., one plate oriented at an arbitrary angle, 0, to the other), the shear strain, y 13 , 
must be set equal to zero for each plate to maintain continuity of rotations. Therefore, if 
Yj 3 is made a fundamental displacement quantity instead of (p,, the shear strain can be set 

equal to zero by simply striking out the appropriate rows and columns in the overall 
stiffness matrix. Performing this operation reduces the stiffness matrix to the same size 
as that for CPT. The VICONOPT code utilizes this procedure for plates that are not 

coplanar. For plates that are coplanar, i.e., 0 = 0, the shear strain in plate 1 is equal to 

that in plate 2. The VICONOPT code handles this situation by creating a substructure 
using the two plates with all DOF present and eliminating the extra DOF before assembly 
into the final stiffness matrix. 

The use of the shear strain, y 13 , as a fundamental displacement quantity requires that 

the effective transverse-shear force per unit length, q 2 , be modified. The modified 
expression for q 2 is obtained from the natural boundary conditions for an edge 
§ 2 =constant, that are derived from the virtual work expression, Eq. (2.8), when y 13 is used 
as a fundamental displacement variable. The modified expression is obtained by 

replacing <p, with the expression w,j — - y 13 in the boundary integral over in Eq. 

R 1 

(2.8). Performing this substitution, integrating by parts, and following the procedure 
outlined in Section 2.6 yields the following modified definition for q 2 : 

x T ( [uj - NT ( [u 2 -z c <p 2 ]) r i ^ 

q 2 = q 2 - N 12 ^w,j — -j - N 22 ^w , 2 - k — +|rni2 — z c n l2j (1 (3- 7 ) 
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The definition for q 2 given in Eq. (3.7) replaces that given in Eq. (2.17c). Note that the 
term m 12 j which appears in the Kirchhoff shear term of CPT is also present for the case 
of SDPT when y 13 , is used as a fundamental displacement quantity. 


3.3 Derivation of the Curved-Plate Stiffness Matrix 


Throughout this section, reference is made to force quantities. Although these 
quantities are forces and moments per unit length, they are designated forces herein for 
convenience. The first step in implementing the present theory into VICONOPT is to 
derive a stiffness matrix that relates the force quantities along the two longitudinal edges, 

^2 = ± ^ , to the displacements along those edges. The desired force and displacement 

quantities are in the direction of the original (undeformed) coordinates. The 
displacement variables are 


1 u, ' 

u 2 


<l>2 

i Y 13 


(3.8) 


where the shear strain, y 13 , has been introduced as a fundamental displacement quantity 


instead of the rotation, <(),. The force variables that correspond to the displacement 


variables given in Eq. (3.8) are 
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(3.9) 
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Note that the effective forces at the boundaries, defined by Eqs. (2.17a) and (2.17b) and 
Eq. (3.7), are being used as the force quantities since, as discussed in Chapter II, they are 
equal to the forces in the direction of the original (undeformed) coordinates. 

The problem may now be reduced to ordinary differential equations in y by assuming 

that the response of the plate in the longitudinal Sj, -direction varies sinusoidally. 
Therefore, if the displacements and forces in the plate are now considered to be functions 
of | 2 , the variables of Eqs. (3.8) and (3.9) may be written as 

(3.10) 

where 



and X is the half- wavelength of the response in the -direction. Since a sinusoidal 

variation in the -direction is assumed, the vector z will involve the amplitudes of the 

displacement and force quantities. The imaginary number, i, has been used in Eqs. (3.8) 
and (3.9) to account for the spatial phase shift that occurs between the perturbation forces 
and displacements which occur at the edges of the plates during buckling for orthotropic 
plates without shear loading and to result in real plate stiffnesses when using the 
exponential expression of Eq. (3.10). 

The next step in the derivation is to express all unknowns in terms of z. A partially 
inverted form of the constitutive relations, Eq. (2.22), is used to express the required 
quantities as functions of the fundamental variables in d and f or terms that may be 
derived from the fundamental variables. The partially inverted constitutive relations are 
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where the linear portion of e n from Eq. (2.2a) is used 


w 

£ 11 = u l,l + — 

K 1 

The variables k u and cp , were defined in Section 2.2 of Chapter II. The constants lv in 

the first portion of Eq. (3.1 1) are calculated from the A, B, and D matrices defined in Eq. 
(2.25). The constants h 77 , h 78 , and h 88 are shear stiffness terms and are calculated using 
the theory presented in |44|. 

Another requirement of the present derivation is to express the relationship between 
q 2 and q 2 without any ^-derivatives. This expression is 


= 


+ ^12 1 W ’l -- 1 Ri CY1J j + N 22 ( < 1 ) 2 _h 78Yl)+[ m 12 -z c n 12]j 

1 “ N 22^88 


(3.12) 


As with the stability equations, only the linear portion of the strain-displacement 
relations are considered in the present derivation 
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Y 12 = u l,2 +u 2,l 

(3.13a) 

W 

e 22 “ u 2,2 + „ 

k 2 

(3.13b) 

u 2 . 

Y 23 = w >2-— “<t>2 

K 2 

(3.13c) 

k 22 = -<$> 2,2 

(3.13d) 

K 12 = ~ (^1,2 + ^ 2 , 1 ) 

(3.13e) 


The expression for k 12 can be re-written after substituting expressions obtained for (J), and 
<j) 2 from Eqs. (2.5a) and (2.5b) and using the linear portion of e ]2 

k 12 = “ Y2 + — j u 2 +2( 1 ) 2 +“- + Y1,2 (3. 13f) 

Using Eqs. (3.11) and (3.12), the strain displacement equations, Eqs. (3.13a) through 
(3.13d) and (3.13f) and the equations, Eqs. (2.16a) through (2.16c) and (2.18a) and 
(2.18b) are written in terms of the elements of z as 

Tz' = Pz or z'=T'' Pz (3.14) 

where a prime denotes differentiation with respect to % 2 . The square matrix T appears in 

the present study as a result of the inclusion of the effects of N 22 and N 12 in the in-plane 
equilibrium equations. This matrix was shown to be the identity matrix when these terms 
were neglected in [2], The presence of off-diagonal terms in this matrix is a fundamental 
difference between the present theory and that presented in [2], 

The elements of z are now assumed to be given by 
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(3.15) 


zj = cj exp 

where (3 is a characteristic root of the differential equation. The number of values of (3 is 

equal to the order of the differential equation system. Substitution of Eq. (3.15) into Eq. 
(3.14) results in the following equation 



(R - (31) c = 0 (3.16a) 

where 

R = bT _1 P (3.16b) 

and I is the identity matrix. The vector c consists of the ^ of Eq. (3.15). The matrix R is 
obtained by premultiplying P by T" 1 . The eigenvalues of the matrix R are the 
characteristic roots of the differential equation. This matrix is not symmetric; however, it 
can be made real by multiplication or division of appropriate rows and columns by the 
imaginary number, i. The elements of the matrices T and P are given in Appendix A for 
both SDPT and CPT. 

For each eigenvalue of R, there exists an eigenvector, c. A matrix C may be defined 
with columns as the eigenvectors, c, the upper half of each column, denoted a, will be 
associated with displacements, and the lower half, denoted b, will be associated with 
forces. The form of C is therefore 


C = 


bi 


«2 

b 2 



(3. 17) 


The next step in the derivation is to write the amplitudes of the displacements and forces 
at the two edges of the plate. Quantities evaluated at are identified with a 
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superscript 1 and quantities evaluated at ^2 = +~ are identified with a superscript 2 as 
follows: 


1 N 

1 = 2 a jk r kexp 

J k=l 

2 N /i B k \ 

j = k 2 ajkrk « p (— ) 



1 N /-ip 

f j = 2 b Jk r k exp — — 
k=l \ 2 

ff - JUjrtexp^) 


(3.18a) 

(3.18b) 

(3.18c) 

(3.18d) 


where the r k are constants determined from the edge values and N is the order of the 
differential equation. Equations (3.18a)-(3.18d) may be written in matrix form as 



Eliminating r from Eqs. (3.19) and (3.20) yields 


(3.19) 

(3.20) 



= K 


d 1 

d 2 


(3.21) 


where K is the stiffness matrix given by 

K=FE‘ (3.22) 

As for the case of CPT, K is real and symmetric for orthotropic plates without in- 
plane shear loading, and it is Hermitian otherwise. Reference [2] presents a discussion of 
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techniques used to ensure that accurate numerical results for K are obtained from Eq. 
(3.22). 


3.4 The Wittrick- Williams Eigenvalue Algorithm 

A brief discussion of the analysis procedure used in VICONOPT is in order. As 
previously mentioned. VICONOPT uses a specialized algorithm for determining any 
natural frequency or buckling load for any given wavelength [29]. The development of 
this algorithm was necessary because the complex stiffnesses defined in the previous 
section are transcendental functions of the load factor and half wavelength of the 
buckling modes of the structure. The eigenvalue problem for determining natural 
frequencies and buckling load factors is therefore transcendental. 

The iterative analysis procedure used in VICONOPT is described in [36]. For this 
procedure, the plate stiffnesses for a given wavelength are evaluated at a series of trial 
values of the eigenvalue being determined until convergence is attained. This eigenvalue 
is either the load factor for buckling or the natural frequency for vibration, and it is 
different than the eigenvalues of the R matrix of Eq. (3.16b). Unless otherwise specified 
by the user, the default initial trial value used in the VICONOPT code is one. For each 
trial value of the eigenvalue considered, the analysis requires the plate stiffnesses as well 
as the number of eigenvalues that lie below the trial value for the entire plate assembly 
assuming the longitudinal edges of each individual to be clamped. A complete 
description of the eigenvalue algorithm is given in [28], Determining the number of 
eigenvalues exceeded by a plate with clamped edges is very difficult except for very 
specialized cases. Therefore, the procedure developed in [28] is used. This procedure 
subdivides each plate into sub-elements with a small enough width such that none of the 
eigenvalues of the sub-elements with clamped edges lie below the trial value. A sub- 
elements is then used as a substructure and is repeatedly doubled until the original plate 
element is obtained. Using a simple procedure at each doubling step [29], the number of 
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eigenvalues that lie below the trial value for the complete plate is returned. This 
procedure is repeated for each plate in the assembly. Using this information and other 
information obtained from the stiffness matrix of the entire assembly, the total number of 
eigenvalues for the entire plate assembly that lie below the trial value is obtained. An 
iterative procedure is then used to refine the trial value until the desired eigenvalue is 
calculated to within the accuracy required. 

One important piece of information required for the analysis procedure described 
herein is the number of subdivisions required for each plate. As seen in Appendix A, all 
of the terms of the R matrix are proportional to the plate width, b. Therefore, all of the 
eigenvalues of R are proportional to b. Furthermore, it is important to note that an 

eigenvalue equal to n corresponds to buckling or vibration with simply supported 
longitudinal edges. By successively halving the value for b until all the real eigenvalues 
of R are less thanrc, a value for the width of the sub-elements for which no eigenvalues 

lie below the eigenvalue for simply supported edges is determined. This width also 
guarantees that no eigenvalues for the sub-elements lie below the eigenvalue for clamped 
edges. This width is that used in the iterative analysis procedure described previously. 
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CHAPTER IV 


NUMERICAL RESULTS 

In this chapter, numerical results are presented using the newly implemented curved- 
plate analysis capability in VICONOPT. Results from several example problems are 
presented to verify the results obtained with this new capability. A convergence study 
using the segmented-plate approach in VICONOPT is performed for an isotropic 
cylindrical shell subjected to uniaxial compression to identify a suitable number of 
segments to be used when comparing results. Results comparing the computational effort 
required by the new analysis to that of the analysis currently in the VICONOPT program 
are also presented for this example. Comparisons of results for several curved plates 
analyzed in Ref. [32] are then made. The effects of including terms related to in-plane 
transverse loads in the in-plane stability equations are examined using a long cylindrical 
tube subjected to in-plane transverse loading. Finally, the curved-plate analysis is used to 
conduct a design-optimization study of a honeycomb-sandwich cylindrical shell 
subjected to uniaxial compression. Comparisons of analyses using both physical and 
tensorial strain measures are made for selected examples, and, where appropriate, results 
using CPT and SDPT are compared. 

4.1 Convergence of the Segmented-Plate Approach 

The convergence of results using the segmented-plate approach in VICONOPT is 
examined for the case of an aluminum cylindrical shell subjected to uniaxial 
compression, see Figure 4.1. The values of the material properties used for this example 

are E = 10.0 x 10 6 psi and v 12 = 0.33. The wall thickness, t, is 0.1 in., and the radius, R 2 is 

60 in. As shown in Reference [49], the critical value for the stress resultant, Nn cr , for 
the axisymmetric buckling of a long isotropic cylindrical shell is 
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(4.1) 


N llcr = 


EC 


R' 


5(l -vf 2 ) 


= 1019.354 lb /in. 


For v 12 = 0.33, the critical half wavelength, k cr , for axisymmetric buckling is shown in 
[49] to be 


^QJ. jt 


l 


2 2 

R2t . = 1.74-y/R^t = 4.255 in. 


12 1-vfa 


(4.2) 


Results illustrating the convergence of the VICONOPT segmented-plate results for 
N]] as a function of the number of segments used to approximate the cylinder are 
shown in Figure 4.2. In this figure, the results of the segmented-plate analysis are shown 
as the solid curve. The theoretical value obtained from Reference [49] is shown as the 
dashed horizontal line. The value obtained using the present curved-plate analysis with 
two curved-plate elements is shown as the open symbol. All results in this figure are 

calculated for the value of X cr given in Eq. (4.2). The VICONOPT results presented in 

this figure are obtained using CPT with tensorial strain measures. As shown in Figure 
4.2, the segmented-plate results converge to the theoretical value when 120 segments are 
used. Therefore, to ensure that converged results are obtained when the segmented-plate 
approach is used to analyze the remaining example problems, sixty segments will be used 
when analyzing curved plates with an included angle of 180 degrees or less, and 120 
segments will be used when analyzing full cylinders. 

This example problem is also used to study the computational requirements of the 
new curved-plate analysis in relation to the segmented-plate approach. A plot of 
normalized CPU time as a function of the number of plate segments used in the 
approximation is shown in Figure 4.3 for the segmented-plate analysis using either CPT 
or SDPT. The normalized CPU time shown in this figure is the CPU time required for 
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the segmented- plate analyses normalized with respect to the CPU time required for the 
curved-plate analysis using two curved-plate elements. The horizontal dashed line is at a 
normalized CPU time of 1.0, and it represents the curved-plate analysis results. As 
shown in the figure, to achieve convergence with 120 flat-plate segments requires 
approximately 3.5 times and 16.7 times as much CPU time as the curved-plate analysis 
for CPT and SDPT, respectively. (For the analysis using SDPT, G 12 = G 13 = G 23 ). One 
consideration to note at this time is that the segmented-plate analysis in VICONOPT is 
implemented to handle the general case of variable geometry, stiffness and loading in the 

^-direction. This approach is therefore not as computationally efficient as it could be for 

the case of constant curvature, stiffness, and loading in that direction (as is the case for 
the curved-plate analysis). One approach to determining the additional computational 
efficiency that may be obtained with the segmented-plate analysis involves defining a 
single, small flat plate that is repeatedly doubled using the substructuring capability in 
VICONOPT until the curved-plate segment is obtained. This technique is referred to 
herein as ‘doubling’. Results relating the computational effort of this approach to the 
curved-plate analysis indicate that further reduction in the computational effort required 
for the segmented-plate analysis can be obtained using this technique. This result occurs 
because the in-plane and out-plane equations are uncoupled in the segmented-plate 
analysis, and analytical expressions for the plate stiffnesses can be used. However, this 
approach is currently not automated in the VICONOPT code, and a separate ‘doubling’ 
effort would have to be made for every curved-plate segment in any given analytical 
model. 

4.2 Buckling of Curved Plates With Widely Varying Curvatures 

The example problems presented in the next two sections are taken from Ref. [32], 
and they are used to verify the results obtained using the new curved-plate analysis in 
VICONOPT. The positive sense of the applied in-plane loadings to be considered in all 
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of the following examples is given in Figure 4.4. The first example problem considered 
studies the buckling of a symmetrically laminated curved plate with clamped longitudinal 
edges as a function of the curvature of the plate. The geometry of the plate is shown in 
Figure 4.5. As shown in the figure, the plate is constructed from symmetrically 
laminated boron/epoxy plies with a [0/90/±45] s layup. To allow for direct comparison of 
results with those presented in [32], the SI units are used for this example and the 
example in the following section. The material properties for a boron/epoxy ply are 
given in Table I. 

The following loadings are considered for this example problem: N n only, N 22 only, 
N 12 only, and combined N„ = N 22 = N 12 . The buckling of this plate subject to these four 
different loadings was investigated while varying the value of the curvature parameter, 


— — , from 1 to 1000. The values of b, R 2 , and (3 used for these analyses are summarized 


in Table 2. Both physical and tensorial strains are used with the new curved-plate 
analysis, while physical strains only are used for the segmented-plate analysis. The 
analysis of [ 32] uses physical strains. All analysis results presented in this section are for 
CPT. The terms involving N 22 and N 12 are included in the in-plane stability equations for 
all analyses. The results of this study are presented in Table 3 for N u loading, in Table 4 
for N 22 loading, in Table 5 for N 12 loading, and in Table 6 for combined N n = N 22 = N 12 
loading. The critical values of the stress resultants presented in these tables were 

X X 

calculated using for the values of — given in the tables. These values of — were 

b b 

presented in Reference [32], The critical values of these stress resultants are also plotted 
as a function of the curvature parameter in Figure 4.6 through Figure 4.9, respectively. 
As shown in these tables, the present analysis compares very well with that presented in 
[32] and with the segmented-plate analysis for widely varying values of the curvature 
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parameter. As shown in Tables 3 - 6, there is no appreciable difference in the results 
obtained using physical and tensorial strains. 

4.3 Buckling of an Unsymmetrically Laminated Curved Plate 

This example problem includes the effect of bending-stretching coupling and shear- 
extension coupling on the buckling of an unsymmetrically laminated curved plate with 
simply supported longitudinal edges. The geometry of the curved plate is shown in 
Figure 4. 10. As shown in the figure, the laminate being studied consists of a 0.0508-cm.- 
thick layer of 2024 aluminum that is reinforced on the inner surface with pairs of ±45’ 
boron/epoxy plies. The material properties for 2024 aluminum are given in Table 1. For 
this example, the number of pairs of ±45° boron/epoxy plies is increased from one to 
seven. The analyses used for this example are identical to those used for the previous 
example. The critical values for buckling of the stress resultants N n and N 12 are 
presented in Tables 7 and 8, respectively. These values are also plotted as a function of 
the number of boron/epoxy plies used in the laminate in Figure 4.11 and Figure 4.12, 
respectively. The agreement between all the analyses is very good. As with the previous 
example, there is no appreciable difference in the results obtained using physical and 
tensorial strains. Results were also computed using SDPT. However, for the case of 
seven pairs of pairs of ±45° boron/epoxy plies, the R/t ratio is still approximately 300, 
and the effects of transverse- shear deformation are minimal. Therefore, as expected, the 
critical values for buckling of the stress resultants N u and N 12 were slightly less than 
those for CPT, but the differences were less than 0.2 percent. With regards to the CPU 
time requirements for this example, the segmented-plate analysis using SDPT required 
approximately 35 times as much CPU time as the curved-plate analysis for the case of 14 
boron/epoxy plies. Furthermore, the results obtained using the 'doubling’ approach 
described in Section 4.1 indicate that the computational efficiencies offered by that 
approach were not realized for this example problem. This result occurs because the 
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coupling that occurs due to the unsymmetric laminate prevents the in-plane and out-of- 
plane equations from being uncoupled, and the same numerical approach for calculating 
the plate stiffnesses as that used for the curved-plate analysis must be used. 

4.4 Effect of N <22 Terms in the In-Plane Stability Equations 

As stated previously, the original segmented-plate analysis in the VICONOPT code 
neglects the effects of the terms involving N 22 and N 12 in the in-plane stability equations. 
This example problem illustrates the effect these terms may have on the buckling of an 
isotropic (aluminum) long cylindrical tube subjected to uniform external pressure. The 
material properties in English units for aluminum are given in Table 9. The geometry of 
this example problem is shown in Figure 4.13. As shown in the figure, only half of the 
tube is modeled since the buckling mode being studied is symmetric (i.e., two full waves 
in the circumferential direction). The pressure load is modeled as an applied N 22 hoop 
loading. The value of the external pressure that would generate this hoop load is obtained 
from the following expression [49] 


P = 


N 22 

r 2 


(4.3) 


Simitses [50] presents a detailed discussion of the buckling of a thin circular ring 
uniformly compressed by external pressure. When considering the behavior of the 
pressure load as the ring buckles, Simitses describes three possible cases. In Case 1 , the 
pressure load is assumed to remain normal to the deflected surface. This loading is 
referred to as a live pressure load. In Case 2, the pressure load is assumed to remain 
parallel to its original direction. This loading is referred to as a dead pressure load. In 
Case 3, the pressure load is assumed to always be directed toward the center of curvature 
of the ring. This loading is referred to as a centrally directed pressure load. Only Cases 1 
and 2 will be discussed in the present report. 
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In Reference [50], Simitses presents in-plane and out-of-plane stability equations for 
the thin circular ring which may be rewritten in the notation used in the present report as 

n 22,2 + ^7" Pcr ( w ’2 “^] "P 2 =0 ( 4 -4) 

02,2 - V 1 - Pcr( R 2 w >22 “ u 2,2 ) " P3 = 0 ( 4 -5) 

k 2 

where p 2 and p 3 are the perturbation values of the applied pressure load in the buckled 
state in the | 2 - and ^-directions, respectively. For the case of a live pressure load in 

which the applied pressure is assumed to remain normal to the deflected surface, q and 
p 3 are (for small deformations) 


P2 = -Pci| w ’2-^j and P3=° ( 4 - 6 ) 

For the case of a dead pressure load in which the applied pressure is assumed to remain 
parallel to its original direction, p 2 and p 3 are 

P 2 = 0 and P 3 = 0 (4.7) 

Substituting Eqs. (4.3) and (4.6) into Eqs. (4.4) and (4.5), yields the following stability 
equations for the case of live pressure loading: 


n 22,2 + ~ = 0 

R 2 


(4.8) 


92,2 


n 22 

r 2 


-N 


22 


cr 



u 2,2 ' 

R 2 / 


= 0 


(4.9) 
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Substituting Eqs. (4.3) and (4.7) into Eqs. (4.4) and (4.5), yields the following stability 
equations for the case of dead pressure loading: 


n 22,2 + 


q2_ 

Rt 


N 


22 cr I w U 2 

W , 2 -■ 


R 


R 


= 0 


(4.10) 


2 / 


02,2 --j^-- N 22 cr 


W,22 - 


u 2,2 

R 


= 0 


(4.11) 


2 / 


For live pressure, the critical value of pressure is shown in [50] to be 


Per = 



Therefore, 


(4.12) 


N 


22 cr 



(4.13) 


For dead pressure, the critical value of pressure is shown in [50] to be 


4EI 



(4.14) 


Therefore, 


N 


22 cr 


4EI 

R2 


(4.15) 
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As discussed in Reference [49], results for the case of a long cylindrical tube uniformly 
compressed by external pressure, may be obtained by considering an elemental ring of 

unit width and using Eqs. (4.12) through (4.15) with E replaced by E/(l-v 2 ) and I 

replaced by t 3 /12. 

External pressure loads are not included in the present analysis. However, an 
equivalent N 22 loading may be calculated using Eq. (4.3). The present analysis treats the 
applied loads as dead loads since no effort is made to modify the applied loads as the 

plate deforms. The stability equations in the ^ 2 and directions, Eqs. (2.11b) and 
(2.1 lc) for the present analysis are written for a thin circular ring subjected to N 22 loading 
by ignoring any terms that involve N I2 , or derivatives with respect to These 
equations are 


n 22,2 


92 N 


22 


R- 


R- 


W, 2 - 


U2_ 

R? 


-FN 


™22 


W , 2 

R? 


+ u 2,22 


= 0 


(4.16) 


92,2 — — N 22 f w ^22 


R' 


R- 


FN 

R- 


22 


w 

r7 


+ u 2,2 


= 0 


(4.17) 


Comparing Eqs. (4.9), (4.11), and (4.17), reveals that if physical strains are used in the 
present analysis (i.e., F = 0), the out-of-plane stability equation is identical to that given 
by Simitses for both live and dead pressure loads. Furthermore, the in-plane stability 
equation for the live pressure load case is recovered by the present analysis if the N 22 
term is neglected in Eq. (4. 16). The dead pressure load case is seen to be recovered when 
the N 22 term is included in Eq. (4.16). Comparing Eqs. (4.6), (4.7), and (4.16), shows that 
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for the case of physical strains, the remaining term involving N 22 is actually equivalent to 
the component of a live pressure load in the § 2 direction (see Eq. (4.6)). 

Buckling results for this example are presented in Table 10. The VICONOPT results 
presented in this table are for physical strains. As previously discussed, the VICONOPT 
result when the N 22 term is neglected in the in-plane stability equation corresponds to the 
case of live pressure load, and the VICONOPT result when the N 22 term is included in the 
in-plane stability equation corresponds to the case of dead pressure load. The results for 
physical strains for the segmented-plate analysis always equal those for the case of dead 
pressure load since the N 22 term in the in-plane stability equation also involves 1/R 2 and it 
therefore drops out of that equation altogether. These results illustrate the dramatic effect 
that the N 22 and N 12 terms in the in-plane stability equations can have on the buckling 
results for curved plates. 

4.5 Design Optimization of a Cylindrical Shell Subject to Uniaxial 

Compression 

The final example utilizes the new curved-plate analysis with the design optimization 
capability of VICONOPT to perform a structural optimization of two different cylindrical 
shell concepts subject to uniform axial compression (N n loading). The two concepts are 
solid-wall construction and honeycomb-sandwich construction. The geometry of this 
example problem is shown in Figure 4. 14. As shown in the figure, the facesheets of the 
honeycomb-sandwich concept are aluminum, and the core is Korex™ aramid paper 
honeycomb core [51]. The solid-wall concept is aluminum. The material properties used 
for the facesheets and core are presented in Table 9. 

Before results for this example are presented, a discussion of the modeling technique 
used to model this cylinder is presented. An analysis of a complete cylinder was 
performed using only one plate element with the new curved-plate analysis capability in 
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VICONOPT. The cylinder is modeled by defining a repetitive cross-section in which 
there is only one node that is connected to itself. However, numerical problems were 
encountered when analyzing closed cylinders with very small wavelengths. The 
following procedure was used to avoid this problem. First, a 45° arc segment is defined. 
Second, a 90° arc segment is defined as a substructure by connecting the original 45° arc 
segment to itself. Similarly, a 180° arc segment is constructed from two 90° arcs. 
Finally, a 360° arc is constructed from two 180° arcs. This substructure is then used to 
define the repetitive cross-section of the cylinder as previously discussed. This modeling 
technique is used for all closed cylinders analyzed in the present report, and no numerical 
problems were encountered when using this technique. 

The design variables for the structural optimization are the thicknesses of the 
facesheets and the core for the sandwich concept and the wall thickness for the solid- wall 
construction. There is no minimum-gage restriction for any design variables. The 
nominal values for these variables are 0.1 in., 0.5 in., and 0.1 in., respectively. The 
design constraints are that the strain in the facesheets or the solid wall cannot exceed 

0.005 in/in and that the stress in the core cannot exceed 1 15 psi in the ^-direction and 55 

psi in the § 2 -direction. The results of this study, including the mass of the optimized 

cylinder and the final values of the design variables are given in Table 11 for the 
honeycomb-sandwich concept and in Table 12 for the solid-wall concept. Results from 
both CPT and SDPT are given in these tables. The optimized mass values are also 
plotted as a function of the applied loading in Figure 4.15. As seen in the tables and the 
figure, the values for optimized mass obtained using CPT are slightly less than those for 
SDPT for the honeycomb-sandwich cylinder as the applied loading is increased. 
However, the values for core thickness obtained using CPT are significantly less than 
those for SDPT for the honeycomb-sandwich cylinder as the applied loading is increased. 
These results are expected because CPT results in an overly stiff approximation since the 
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effects of transverse-shear deformation are neglected. This overly stiff approximation 
results in higher buckling loads for a given core thickness. Therefore, the core thickness 
and the optimum mass obtained is less than that obtained using SDPT. The optimized 
mass values for the solid-wall construction are much greater than those for the 
honeycomb-sandwich construction. The results for CPT and SDPT are nearly identical 
for the solid-wall construction, as expected. 
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CHAPTER V 


CONCLUDING REMARKS 

The VICONOPT computer code is an exact analysis and optimum design program 
that includes the buckling and vibration analyses of prismatic assemblies of flat, in-plane- 
loaded anisotropic plates. In the present report, the capability to analyze curved-plate 
segments exactly has been added to the VICONOPT code. Non-linear curved-plate 
equilibrium equations have been formulated using the principle of virtual work, and 
linearized stability equations that describe the response of the plate just after buckling 
occurs were derived following the application of several simplifying assumptions. 
Finally, modifications to these equations were made to allow the reference surface of the 
plate to be located at a distance z c from the centroidal surface. 

The analysis methodology described in the present report improves upon the existing 
methodology in the VICONOPT code (referred to herein as the segmented-plate analysis) 
which requires that curved-plate segments be subdivided into several flat-plate elements 
that must be subsequently joined at their longitudinal edges to approximate the curved- 
plate geometry. The new analysis formulation allows either classical plate theory (CPT) 
or first-order shear deformation plate theory (SDPT) to be used. Furthermore, anisotropic 
laminates having fully populated A, B, and D stiffness matrices may be analyzed. The 
analysis described in the present report is an example of an exact finite-strip method 
(FSM) since it uses a stiffness matrix that is derived by direct solution to the stability 
equations. 

One additional capability that has been incorporated into the VICONOPT code as part 
of the present report is the option to use plate elements (flat or curved) that are based 
upon nonlinear strain-displacement relations that contain terms from either physical or 
tensorial strain measures. A second capability that has been added is the ability to 
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include the effect of terms associated with in-plane transverse and in-plane shear loading 
in the in-plane stability equations. The original VICONOPT code neglects these terms. 

Results from the present curved-plate analysis capability compare very well with a 
closed-form solution and the existing segmented-plate analysis for the buckling of a long 
isotropic cylinder. The present analysis also compares well with results from the 
literature for symmetrically laminated curved plates with widely varying curvatures and 
with unsymmetrically laminated plates that include the effect of extensional-bending and 
shear-extension coupling. No appreciable effects of using tensorial versus physical 
strains are noted in these examples. The present curved-plate analysis was also shown to 
require significantly less computational effort than the segmented-plate analysis. An 
alternate approach for the segmented-plate analysis that offers additional computational 
savings for certain classes of problems has been investigated. However, this approach 
requires greater user effort, and it is currently not implemented in the VICONOPT code. 

A significant effect of either including or neglecting the terms associated with an 
applied in-plane transverse loading (i.e., N 22 loading) in the in-plane stability equations 
was noted when analyzing a long cylindrical tube subjected to uniform external pressure. 
The symmetry of the buckling mode for this problem allowed it to be modeled as a half 
cylinder, and the pressure load was simulated with an equivalent hoop (N 22 ) loading. The 
buckling results for this problem were shown to change by a factor of 3/4 when the terms 
associated with the N 22 loading were neglected in the in-plane stability equations. This 
result illustrates the effect that the treatment of the in-plane stability equations can have 
on the buckling results for curved plates. 

Finally, the present curved-plate analysis was used to conduct a design-optimization 
study of two different cylindrical shells subject to uniform axial compression (N n 
loading). One shell was constructed from a honeycomb-sandwich wall construction, and 
the other was a solid-wall construction. The values of mass for the optimized solid-wall 
design were consistently higher than those for the honeycomb-sandwich construction. 
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However, there was no difference between results using CPT and SDPT for the solid-wall 
cylinder. The values of mass for the optimized honeycomb-sandwich cylinder using CPT 
were slightly less than those for SDPT as the applied loading was increased. However, 
the values of core thickness for the optimized honeycomb-sandwich cylinder using CPT 
were significantly less than those for SDPT as the applied loading was increased. This 
trend occurred because CPT results in an overly stiff approximation since the effects of 
transverse-shear flexibility are neglected. This overly stiff approximation results in 
higher buckling loads and, thus, a lower optimum mass. 

One area for future work includes retaining the curvature terms in the longitudinal 
direction and implementing the capability to analyze shells of revolution. The analysis 
can also be modified to allow vibration analyses to be performed. Another enhancement 
that can be made to the present analysis is to remove the restriction that when the terms 
associated with in-plane transverse and in-plane shear loading are retained in the in-plane 
stability equations, the centroidal surface and reference surface must coincide. 
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APPENDIX A 


MATRICES FOR DETERMINING CHARACTERISTIC ROOTS 

The eigenvalues of matrix R in Eq. (3.16a) are the characteristic roots of the 
differential equations describing the behavior of the plate. The 10-by-10 matrix R is 
calculated from the matrices T and P as shown in Eq. (3.16b). The non-zero elements of 
the T matrix are 
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The underlined terms given above and subsequently are those terms that drop out of the 
equation when the effects of N 22 and N 12 in the in-plane stability equations are neglected. 
The 11011 -zero elements of the P matrix are 
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s = 1 - h88 N 22 

The expressions for the elements of the T and P matrices for the case of classical 
plate theory are obtained by setting the transverse-shear strains, y 13 and y 23 , equal to zero 

and using the resulting expression <t>j = w j (recall that — equals to zero). The 

R 1 

partially inverted stress-strain relations given in Eq. (3.11), are modified such that m 12 
and k 12 are interchanged. For the classical case, only four stability equations, Eqs. 

(2.16a), (2.16b), (2.16c), and (2.18a) are used since Eq. (2.18b) is satisfied by 
incorporation into the final form of Eq. (2.16c). The same steps used for the transverse- 
shear case are followed to generate T and P matrices of order eight. The elements of the 
T and P matrices with a superscript * given previously for the transverse-shear case also 
apply for the classical case if 1 is subtracted from any index greater than 4. The non-zero 
elements of the T matrix that are not given in the results for transverse shear are 
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The non-zero elements of the P matrix that are not given in the results for transverse 
shear are 
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Table 1. Material properties for boron/epoxy plies and 2024 aluminum (SI units). 


Material 

EnxlO- 10 , 

N/m 2 

E22X10-J0, 

N/m 2 

Gi 2 xlO-10, 

N/m 2 

V12 

p, kg/m 3 

Boron/epoxy 

20.69 

1.86 

0.48 

0.21 

2006.8 

Aluminum 2024 

7.38 

7.38 

2.76 

0.33 

2768.0 


Table 2. Geometric parameters used to vary the curvature parameter, . 

R2t 


b 2 

R 2 t 

b, cm. 

R 2 , cm. 

(3, degrees 

1 

24.4002 

5760.3570 

0.25264 

5 

24.4005 

1152.1180 

1 .2632 

10 

25.4020 

576.12905 

2.5262 

30 

25.4185 

192.2917 

7.5738 

50 

25.4513 

115.67310 

12.6067 

100 

25.6036 

58.53098 

25.0633 

300 

27.1026 

21.86161 

71.0315 

500 

29.6186 

15.66554 

108.3281 

700 

32.6900 

13.63059 

137.4115 


37.7873 

12.75046 

169.8018 


71 

























































C3 



72 



Table 4. Critical value of stress resultant N 22 for buckling of a symmetrically laminated long curved plate with clamped longitudinal 
edges. 

VICONOPT a , N 22cr ,N/m 

Physical strains Tensorial strains 

Curved-plate theory Segmented-plate theory Curved-plate theory 

5 856 

8 077 

11 780 

11 767 

11 743 

11 632 

10 607 

9 143 

7 696 

5 887 

' N 22 terms included in the in-plane stability equations. 
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Table 9. Material properties for aluminum and Korex™ honeycomb core (English Engineering units). 

EnxlO' 6 , E22XIO' 6 , Gi2xl0~ 6 , Gi3xl0~ 6 , G23XIO' 6 , V12 p, lb/in 3 


No 




a 

3 


_a 

23 


.3 


.a 

23 


3 
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Table 12. Design-optimization results for a solid- wall cylinder subjected to Nn loading. 

Transverse shear plate theory (physical strains) 
t wa ib in. mass, lb m 

924.69 

1296.8 

1622.0 

1873.9 

2082.2 

2931.9 

3697.7 

9\6Zt 

04 

O 

o 

cn 

© 

0.179 

0.207 

0.230 

0.324 

0.409 

0.474 

Classical plate theory (physical strains) 
t wa ib in- mass, lb m 

924.68 

1296.8 

1622.0 

1873.8 

2082.1 

2931.8 

3697.5 

4291.3 
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Figure 1.1 Typical longitudinally stiffened plate structures. 




Figure 1.2 Segmented representation of curved-plate geometry currently used by 
VICONOPT. 
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Figure 2.2 Sign convention for applied in-plane loads and relation of reference surface 
to centroidal surface. 
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Plate 1 


Plate 2 



(a) plates 1 and 2 coplanar 



(c) plate 2 rotated to arbitray angle, 0 
Figure 3.1 Displacements and rotations at a typical plate junction. 
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Segmented-plate theory CPU times normalized 
with respect to curved-plate theory time 


Normalized 
CPU time 



Figure 4.3 Normalized CPU time requirements for the segmented-plate approach as a 
function of the number of segments used in the approximation. 
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Figure 4.6 Critical value of stress resultant N n for buckling of a symmetrically 
laminated curved plate with clamped longitudinal edges. 




. ^ 


Nl2 



Figure 4.8 Critical value of stress resultant N 12 for buckling of a symmetrically 
laminated curved plate with clamped longitudinal edges. 


90 







Curved-plate theory, physical strains 
Segmented-plate theory, physical strains 
Curved-plate theory, tensorial strains 
Reference [32] 


2 34567 


2 34567 


2 3 4 5 6 7 


Curvature parameter, b 2 /R 2 t 


Figure 4.9 Critical value of stress resultant N n =N 22 =N 12 for buckling of a symmetrically 
laminated curved plate with clamped longitudinal edges. 








Figure 4.10 


ui and w = 0 



Section A-A 


+45° B/E 
-45° B/E 
+45° B/E 
-45° B/E 



* 0.014 cm. 

try (typical) 


0.0508 cm. 


Unsymmetrically laminated aluminum and boron/epoxy (B/E) curved 
plate with simply supported edges subjected to applied in-plane loads. 
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Number of boron/epoxy plies 

Figure 4. 11 Critical value of stress resultant N n for buckling of an unsymmetrically 

laminated aluminum and boron/epoxy (B/E) curved plate with simply 
supported longitudinal edges. 
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Nl2 



Figure 4.12 Critical value of stress resultant N 12 for buckling of an unsymmetrically 
laminated aluminum and boron/epoxy (B/E) curved plate with simply 
supported longitudinal edges. 
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Section A-A 
Solid wall 

1 0. 1 in. nom. 

=^=T 


Aluminum 


Section A-A 



facesheets honeycomb core 


Figure 4.14 Cylindrical shell subjected to uniform axial compression (N n loading). 
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Figure 4.15 


Optimized cylinder mass as a function of the applied loading for a 
cylindrical shell. 
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